DifferentialEquations.jl Workshop Exercise Solutions
Chris Rackauckas
using DifferentialEquations using Sundials using BenchmarkTools using Plots
Problem 1: Investigating Sources of Randomness and Uncertainty in a Biological System
Part 1: Simulating the Oregonator ODE model
using DifferentialEquations, Plots function orego(du,u,p,t) s,q,w = p y1,y2,y3 = u du[1] = s*(y2+y1*(1-q*y1-y2)) du[2] = (y3-(1+y1)*y2)/s du[3] = w*(y1-y3) end p = [77.27,8.375e-6,0.161] prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,360.0),p) sol = solve(prob) plot(sol)
plot(sol,vars=(1,2,3))
Part 2: Investigating Stiffness
using BenchmarkTools prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,50.0),p) @btime sol = solve(prob,Tsit5())
1.217 s (8723143 allocations: 920.67 MiB)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 872306-element Array{Float64,1}:
0.0
0.01618926718934831
0.02355386004837834
0.03818038870154586
0.050503515877727514
0.06810672932191658
0.08676359998206734
0.11145368602241688
0.14105967462147356
0.18104879156165962
⋮
49.99977330536325
49.99980456142745
49.999835817515255
49.999867073624586
49.999898329755446
49.99992958590576
49.99996084207554
49.999992098266844
50.0
u: 872306-element Array{Array{Float64,1},1}:
[1.0, 2.0, 3.0]
[1.7128564042197614, 1.9996098373795999, 2.9959141611121862]
[1.8376268914687968, 1.9993653073090198, 2.994474646468457]
[1.9480445809808178, 1.9988333244430836, 2.991907642632475]
[1.9807789479174538, 1.998364632682339, 2.989876120098015]
[1.996520358969301, 1.9976843022063284, 2.9870473687154533]
[2.0012471416469095, 1.9969587120867922, 2.9840850652644586]
[2.003267094253373, 1.9959962346456372, 2.980190667568818]
[2.0046071951018165, 1.9948405279663373, 2.9755485736940304]
[2.0062040975915965, 1.9932773146432707, 2.969322732597494]
⋮
[1.00114451241949, 1453.0173573419604, 414.83224206133156]
[1.0011445128905938, 1453.0163492345089, 414.8301595725294]
[1.001144513536549, 1453.0153411262695, 414.82807709263454]
[1.001144514166616, 1453.014333017309, 414.82599462178484]
[1.0011445147807905, 1453.013324907627, 414.8239121599803]
[1.0011445151883325, 1453.0123167972909, 414.82182970735875]
[1.0011445153892404, 1453.0113086863003, 414.8197472639202]
[1.001144515574252, 1453.0103005745884, 414.8176648295267]
[1.0008765717435082, 1453.0100456736175, 414.8171383809634]
@btime sol = solve(prob,Rodas5())
478.942 μs (1909 allocations: 130.16 KiB)
retcode: Success
Interpolation: 3rd order Hermite
t: 110-element Array{Float64,1}:
0.0
0.019615259849088615
0.029598314714131158
0.04705295553350644
0.06489958093933189
0.08933251171067431
0.12069400166576917
0.16655311655246774
0.24089140897016648
0.39558909491704786
⋮
26.756905610888992
27.982111658219903
29.768997154114096
32.21837697976615
35.093850201346655
38.49798110093118
42.33811919585127
46.60842194880463
50.0
u: 110-element Array{Array{Float64,1},1}:
[1.0, 2.0, 3.0]
[1.7804115041903392, 1.9994992840408727, 2.995224421497252]
[1.898773632635922, 1.9991507098568697, 2.9933805881501456]
[1.9745775749460168, 1.9984968888022705, 2.9904382700551317]
[1.9949959346655894, 1.9978087397951183, 2.9875591897227847]
[2.0015958931121642, 1.9968586608477479, 2.98367866834122]
[2.003748190575679, 1.9956356930368464, 2.9787387129953866]
[2.0056429388535917, 1.9938442509772465, 2.9715736894920433]
[2.0085949421229565, 1.9909335157971781, 2.960099467726684]
[2.014815188384092, 1.9848502001186519, 2.936770263171178]
⋮
[1.0009510454262696, 1052.1681949981978, 17454.97704553619]
[1.000790082105047, 1266.4223517298105, 14330.342720311946]
[1.0006713873660182, 1490.2781714142227, 10747.93771088393]
[1.000598803847115, 1670.9447027478102, 7245.705166049239]
[1.000568993307521, 1758.4723173221284, 4560.988616721742]
[1.000569273504183, 1757.6100577789323, 2636.982996349979]
[1.000594225030407, 1683.8471494545056, 1421.4818618119598]
[1.0006409946157637, 1561.0560213127278, 715.2527024515273]
[1.0006887475677544, 1452.8969192375328, 414.7220773988324]
(Optional) Part 3: Specifying Analytical Jacobians (I)
(Optional) Part 4: Automatic Symbolicification and Analytical Jacobian Calculations
Part 5: Adding stochasticity with stochastic differential equations
function orego(du,u,p,t) s,q,w = p y1,y2,y3 = u du[1] = s*(y2+y1*(1-q*y1-y2)) du[2] = (y3-(1+y1)*y2)/s du[3] = w*(y1-y3) end function g(du,u,p,t) du[1] = 0.1u[1] du[2] = 0.1u[2] du[3] = 0.1u[3] end p = [77.27,8.375e-6,0.161] prob = SDEProblem(orego,g,[1.0,2.0,3.0],(0.0,30.0),p) sol = solve(prob,SOSRI()) plot(sol)
sol = solve(prob,ImplicitRKMil()); plot(sol)
sol = solve(prob,ImplicitRKMil()); plot(sol)
Part 6: Gillespie jump models of discrete stochasticity
Part 7: Probabilistic Programming / Bayesian Parameter Estimation with DiffEqBayes.jl + Turing.jl (I)
The data was generated with:
function orego(du,u,p,t) s,q,w = p y1,y2,y3 = u du[1] = s*(y2+y1*(1-q*y1-y2)) du[2] = (y3-(1+y1)*y2)/s du[3] = w*(y1-y3) end p = [60.0,1e-5,0.2] prob = ODEProblem(orego,[1.0,2.0,3.0],(0.0,30.0),p) sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
retcode: Success
Interpolation: 3rd order Hermite
t: 48825-element Array{Float64,1}:
0.0
0.0001377354452002734
0.0002010718419122773
0.0003021998626318717
0.0004033278833514661
0.0005062369821758017
0.0006097468311490204
0.0007142127440433809
0.0008192699078400218
0.0009249477742217067
⋮
29.8029621721172
29.830644478725933
29.858326785334665
29.886009659186307
29.91369253303795
29.94137540688959
29.96906121150096
29.99674701611233
30.0
u: 48825-element Array{Array{Float64,1},1}:
[1.0, 2.0, 3.0]
[1.0082299897608653, 1.9999976854327, 2.9999450200919755]
[1.0119917072902584, 1.999996608412124, 2.999919814481464]
[1.0179684138472678, 1.9999948722835914, 2.99987966834307]
[1.0239089650942212, 1.9999931160074038, 2.9998396435319274]
[1.029917251895588, 1.9999913082422207, 2.999799037595289]
[1.0359233048539687, 1.999989469133448, 2.9997583198166615]
[1.0419471308701034, 1.9999875920280812, 2.9997173524755154]
[1.0479670917315083, 1.9999856831456868, 2.9996762806603012]
[1.0539844577981918, 1.9999837417184405, 2.99963509426064]
⋮
[1.000649179056226, 1541.3915652677524, 2704.0497424834743]
[1.0006492543168537, 1541.2130026445143, 2689.1257668560265]
[1.0006493324193109, 1541.027741407077, 2674.2841889364954]
[1.0006494133487425, 1540.8358216266624, 2659.5242521879404]
[1.0006494970871955, 1540.6372908555195, 2644.845809127133]
[1.0006495836184481, 1540.4321926696598, 2630.248409804933]
[1.0006496729360377, 1540.220547649861, 2615.7300741273025]
[1.0006497650150368, 1540.0024202595741, 2601.291906654189]
[1.00064977601455, 1539.9763674790563, 2599.600715743143]
(Optional) Part 8: Using DiffEqBiological's Reaction Network DSL
Problem 2: Fitting Hybrid Delay Pharmacokinetic Models with Automated Responses (B)
Part 1: Defining an ODE with Predetermined Doses
function onecompartment(du,u,p,t) Ka,Ke = p du[1] = -Ka*u[1] du[2] = Ka*u[1] - Ke*u[2] end p = (Ka=2.268,Ke=0.07398) prob = ODEProblem(onecompartment,[100.0,0.0],(0.0,90.0),p) tstops = [24,48,72] condition(u,t,integrator) = t ∈ tstops affect!(integrator) = (integrator.u[1] += 100) cb = DiscreteCallback(condition,affect!) sol = solve(prob,Tsit5(),callback=cb,tstops=tstops) plot(sol)
Part 2: Adding Delays
function onecompartment_delay(du,u,h,p,t) Ka,Ke,τ = p delayed_depot = h(p,t-τ)[1] du[1] = -Ka*u[1] du[2] = Ka*delayed_depot - Ke*u[2] end p = (Ka=2.268,Ke=0.07398,τ=6.0) h(p,t) = [0.0,0.0] prob = DDEProblem(onecompartment_delay,[100.0,0.0],h,(0.0,90.0),p) tstops = [24,48,72] condition(u,t,integrator) = t ∈ tstops affect!(integrator) = (integrator.u[1] += 100) cb = DiscreteCallback(condition,affect!) sol = solve(prob,MethodOfSteps(Rosenbrock23()),callback=cb,tstops=tstops) plot(sol)
Part 3: Automatic Differentiation (AD) for Optimization (I)
Part 4: Fitting Known Quantities with DiffEqParamEstim.jl + Optim.jl
The data was generated with
p = (Ka = 0.5, Ke = 0.1, τ = 4.0)
(Ka = 0.5, Ke = 0.1, τ = 4.0)
Part 5: Implementing Control-Based Logic with ContinuousCallbacks (I)
Part 6: Global Sensitivity Analysis with the Morris and Sobol Methods
Problem 3: Differential-Algebraic Equation Modeling of a Double Pendulum (B)
Part 1: Simple Introduction to DAEs: Mass-Matrix Robertson Equations
function f(du, u, p, t) du[1] = -p[1]*u[1] + p[2]*u[2]*u[3] du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2] du[3] = u[1] + u[2] + u[3] - 1. end M = [1 0 0; 0 1 0; 0 0 0.] p = [0.04, 10^4, 3e7] u0 = [1.,0.,0.] tspan = (0., 1e6) prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p) sol = solve(prob, Rodas5()) plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))
Part 2: Solving the Implicit Robertson Equations with IDA
# Robertson Equation DAE Implicit form function h(out, du, u, p, t) out[1] = -p[1]*u[1] + p[2]*u[2]*u[3] - du[1] out[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2] - du[2] out[3] = u[1] + u[2] + u[3] - 1. end p = [0.04, 10^4, 3e7] du0 = [-0.04, 0.04, 0.0] u0 = [1.,0.,0.] tspan = (0., 1e6) differential_vars = [true, true, false] prob = DAEProblem(h, du0, u0, tspan, p, differential_vars = differential_vars) sol = solve(prob, IDA()) plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))
Part 3: Manual Index Reduction of the Single Pendulum
Consider the equation:
\[ x^2 + y^2 = L \]
Differentiating once with respect to time:
\[ 2x\dot{x} + 2y\dot{y} = 0 \]
A second time:
\[ \begin{align} {\dot{x}}^2 + x\ddot{x} + {\dot{y}}^2 + y\ddot{y} &= 0 \\ u^2 + v^2 + x(\frac{x}{mL}T) + y(\frac{y}{mL}T - g) &= 0 \\ u^2 + v^2 + \frac{x^2 + y^2}{mL}T - yg &= 0 \\ u^2 + v^2 + \frac{T}{m} - yg &= 0 \end{align} \]
Our final set of equations is hence
\[ \begin{align} \ddot{x} &= \frac{x}{mL}T \\ \ddot{y} &= \frac{y}{mL}T - g \\ \dot{x} &= u \\ \dot{y} &= v \\ u^2 + v^2 -yg + \frac{T}{m} &= 0 \end{align} \]
We finally obtain $T$ into the third equation. This required two differentiations with respect to time, and so our system of equations went from index 3 to index 1. Now our solver can handle the index 1 system.
Part 4: Single Pendulum Solution with IDA
function f(out, da, a, p, t) (L, m, g) = p u, v, x, y, T = a du, dv, dx, dy, dT = da out[1] = x*T/(m*L) - du out[2] = y*T/(m*L) - g - dv out[3] = u - dx out[4] = v - dy out[5] = u^2 + v^2 - y*g + T/m nothing end # Release pendulum from top right u0 = zeros(5) u0[3] = 1.0 du0 = zeros(5) du0[2] = 9.81 p = [1,1,9.8] tspan = (0.,100.) differential_vars = [true, true, true, true, false] prob = DAEProblem(f, du0, u0, tspan, p, differential_vars = differential_vars) sol = solve(prob, IDA()) plot(sol, vars=(3,4))
Part 5: Solving the Double Penulum DAE System
For the double pendulum: The equations for the second ball are the same as the single pendulum case. That is, the equations for the second ball are:
\[ \begin{align} \ddot{x_2} &= \frac{x_2}{m_2L_2}T_2 \\ \ddot{y_2} &= \frac{y_2}{m_2L_2}T_2 - g \\ \dot{x_2} &= u \\ \dot{y_2} &= v \\ u_2^2 + v_2^2 -y_2g + \frac{T_2}{m_2} &= 0 \end{align} \]
For the first ball, consider x_1^2 + y_1^2 = L $
\[ \begin{align} x_1^2 + x_2^2 &= L \\ 2x_1\dot{x_1} + 2y_1\dot{y_1} &= 0 \\ \dot{x_1}^2 + \dot{y_1}^2 + x_1(\frac{x_1}{m_1L_1}T_1 - \frac{x_2}{m_1L_2}T_2) + y_1(\frac{y_1}{m_1L_1}T_1 - g - \frac{y_2}{m_1L_2}T_2) &= 0 \\ u_1^2 + v_1^2 + \frac{T_1}{m_1} - \frac{x_1x_2 + y_1y_2}{m_1L_2}T_2 &= 0 \end{align} \]
So the final equations are:
\[ \begin{align} \dot{u_2} &= x_2*T_2/(m_2*L_2) \dot{v_2} &= y_2*T_2/(m_2*L_2) - g \dot{x_2} &= u_2 \dot{y_2} &= v_2 u_2^2 + v_2^2 -y_2*g + \frac{T_2}{m_2} &= 0 \dot{u_1} &= x_1*T_1/(m_1*L_1) - x_2*T_2/(m_2*L_2) \dot{v_1} &= y_1*T_1/(m_1*L_1) - g - y_2*T_2/(m_2*L_2) \dot{x_1} &= u_1 \dot{y_1} &= v_1 u_1^2 + v_1^2 + \frac{T_1}{m_1} + \frac{-x_1*x_2 - y_1*y_2}{m_1L_2}T_2 - y_1g &= 0 \end{align} \]
function f(out, da, a, p, t) L1, m1, L2, m2, g = p u1, v1, x1, y1, T1, u2, v2, x2, y2, T2 = a du1, dv1, dx1, dy1, dT1, du2, dv2, dx2, dy2, dT2 = da out[1] = x2*T2/(m2*L2) - du2 out[2] = y2*T2/(m2*L2) - g - dv2 out[3] = u2 - dx2 out[4] = v2 - dy2 out[5] = u2^2 + v2^2 -y2*g + T2/m2 out[6] = x1*T1/(m1*L1) - x2*T2/(m2*L2) - du1 out[7] = y1*T1/(m1*L1) - g - y2*T2/(m2*L2) - dv1 out[8] = u1 - dx1 out[9] = v1 - dy1 out[10] = u1^2 + v1^2 + T1/m1 + (-x1*x2 - y1*y2)/(m1*L2)*T2 - y1*g nothing end # Release pendulum from top right u0 = zeros(10) u0[3] = 1.0 u0[8] = 1.0 du0 = zeros(10) du0[2] = 9.8 du0[7] = 9.8 p = [1,1,1,1,9.8] tspan = (0.,100.) differential_vars = [true, true, true, true, false, true, true, true, true, false] prob = DAEProblem(f, du0, u0, tspan, p, differential_vars = differential_vars) sol = solve(prob, IDA()) plot(sol, vars=(3,4)) plot(sol, vars=(8,9))
Problem 4: Performance Optimizing and Parallelizing Semilinear PDE Solvers (I)
Part 1: Implementing the BRUSS PDE System as ODEs
using DifferentialEquations, Sundials, Plots # initial condition function init_brusselator_2d(xyd) N = length(xyd) u = zeros(N, N, 2) for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] u[I,1] = 22*(y*(1-y))^(3/2) u[I,2] = 27*(x*(1-x))^(3/2) end u end N = 32 xyd_brusselator = range(0,stop=1,length=N) u0 = vec(init_brusselator_2d(xyd_brusselator)) tspan = (0, 22.) p = (3.4, 1., 10., xyd_brusselator) brusselator_f(x, y, t) = ifelse((((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) && (t >= 1.1), 5., 0.) using LinearAlgebra, SparseArrays du = ones(N-1) D2 = spdiagm(-1 => du, 0=>fill(-2.0, N), 1 => du) D2[1, N] = D2[N, 1] = 1 D2 = 1/step(xyd_brusselator)^2*D2 tmp = Matrix{Float64}(undef, N, N) function brusselator_2d_op(du, u, (D2, tmp, p), t) A, B, α, xyd = p dx = step(xyd) N = length(xyd) α = α/dx^2 du = reshape(du, N, N, 2) u = reshape(u, N, N, 2) @views for i in axes(u, 3) ui = u[:, :, i] dui = du[:, :, i] mul!(tmp, D2, ui) mul!(dui, ui, D2') dui .+= tmp end @inbounds begin for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] i = I[1] j = I[2] du[i,j,1] = α*du[i,j,1] + B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t) du[i,j,2] = α*du[i,j,2] + A*u[i,j,1] - u[i,j,1]^2*u[i,j,2] end end nothing end prob1 = ODEProblem(brusselator_2d_op, u0, tspan, (D2, tmp, p)) sol1 = @time solve(prob1, TRBDF2(autodiff=false));
9.387620 seconds (9.52 M allocations: 585.980 MiB, 1.31% gc time)
retcode: Success
Interpolation: 3rd order Hermite
t: 79-element Array{Float64,1}:
0.0
1.0522804268928948e-10
1.1575084695821841e-9
1.1680312738511129e-8
1.1690835542780058e-7
1.169188782320695e-6
1.1691993051249639e-5
3.96415424995973e-5
0.0001251040827040069
0.00037992550155953065
⋮
18.252302769116575
18.556761040654674
19.073892463800625
19.428284008927427
20.09879123998184
20.480538570171507
20.999444301927937
21.403330440737882
22.0
u: 79-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 2.5250877095783344,
2.2620667554742258, 1.9735248771761977, 1.665005111992191, 1.34336401668221
03, 1.0172186542655526, 0.6977464117458191, 0.4003323380813969, 0.148922584
53196738, 0.0]
[0.00011790503536898665, 0.00011790503536909906, 0.00011790503536928888, 0
.00011790503536951345, 0.0001179050353697547, 0.00011790503537000097, 0.000
11790503537024383, 0.00011790503537047684, 0.00011790503537069474, 0.000117
90503537089335 … 2.5250585327638797, 2.262041954652914, 1.973505463224665
6, 1.6649923614746842, 1.3433596406307873, 1.0172251412113864, 0.6977678519
884, 0.4003770605722655, 0.14902217713448462, 0.000144701505282153]
[0.0012949971428852717, 0.0012949971428992155, 0.001294997142922781, 0.001
2949971429507092, 0.0012949971429807212, 0.0012949971430113658, 0.001294997
1430415896, 0.0012949971430705825, 0.001294997143097699, 0.0012949971431224
178 … 2.52476680944311, 2.2617940004952306, 1.973311391982252, 1.66486494
7832933, 1.343316013313132, 1.0172902298879962, 0.6979827036906551, 0.40082
59597669359, 0.15001758272059923, 0.0015893132679671673]
[0.012876383252354838, 0.012876383254089876, 0.012876383257051385, 0.01287
6383260614965, 0.012876383264461243, 0.012876383268395485, 0.01287638327227
9336, 0.012876383276007026, 0.01287638327949474, 0.012876383282674856 … 2
.5218540714350746, 2.2593198858887353, 1.971377545998155, 1.663600045556335
, 1.3428932614037994, 1.0179637350768072, 0.7001796234851249, 0.40546907561
279366, 0.159920626066569, 0.015802820975871525]
[0.11431809668801969, 0.11431809701540827, 0.11431809761055266, 0.11431809
838668756, 0.11431809926467854, 0.1143181001852524, 0.11431810110619822, 0.
1143181019968641, 0.1143181028341037, 0.11431810359985839 … 2.49319122864
72744, 2.2351481567086453, 1.9527784378871984, 1.6519892724152478, 1.340282
573226162, 1.0275217594291501, 0.7274071832384641, 0.4600255371073474, 0.25
403236055497325, 0.14029945313886596]
[0.6824546276302552, 0.6824546576303283, 0.6824547160548516, 0.68245479982
66596, 0.6824549045779489, 0.6824550250157978, 0.6824551553829176, 0.682455
2898965248, 0.6824554230873283, 0.6824555500182544 … 2.2673901444639477,
2.0670335967015725, 1.856281803467327, 1.6433610322349588, 1.43773260755362
28, 1.249622628038894, 1.089014709165961, 0.964086299638045, 0.879715279699
5076, 0.8375618491417234]
[1.6813894581935167, 1.6813894979141795, 1.6813895760496735, 1.68138969005
22025, 1.6813898362183566, 1.6813900097580707, 1.681390204842709, 1.6813904
14668493, 1.6813906315812168, 1.681390847292351 … 1.8814741145279639, 1.9
062408314582533, 1.932556304652824, 1.9593615779722884, 1.9854302030584638,
2.009456254944375, 2.030184671316311, 2.0465572620802512, 2.05782082872771
7, 2.0635421902068396]
[1.5465563234604633, 1.5465563151032578, 1.546556298667805, 1.546556274705
925, 1.5465562440308058, 1.5465562077004953, 1.5465561669954782, 1.54655612
33881666, 1.5465560785018047, 1.546556034057256 … 1.9340770272586247, 1.9
289445685108866, 1.9236144878920916, 1.918294922809417, 1.9132050182112592,
1.90856340239867, 1.9045749934046978, 1.9014179583244304, 1.89923295644000
23, 1.8981161877601427]
[1.5715361948410418, 1.5715361953311315, 1.57153619629489, 1.5715361976997
826, 1.5715361994977792, 1.571536201626338, 1.571536204009821, 1.5715362065
614462, 1.5715362091858556, 1.5715362117823382 … 1.9247785400997206, 1.92
53926733129951, 1.926022344707277, 1.926643291184795, 1.9272311943883393, 1
.9277626997857489, 1.9282164458524176, 1.9285740360607109, 1.92882088655398
73, 1.9289468947014914]
[1.5684925853534752, 1.5684925853540483, 1.5684925853552447, 1.56849258535
7125, 1.5684925853597242, 1.5684925853630174, 1.5684925853669094, 1.5684925
853712433, 1.5684925853758194, 1.5684925853804146 … 1.9259264057842067, 1
.9258966605641663, 1.9258662870639163, 1.9258364529304535, 1.92580831012103
66, 1.925782949722734, 1.9257613584940534, 1.9257443790768893, 1.9257326756
936477, 1.925726706942093]
⋮
[3.519620884001627, 3.5196209320997647, 3.5196209820910895, 3.519621032110
1122, 3.5196210801507464, 3.5196211241431548, 3.5196211620586335, 3.5196211
92036574, 3.5196212125280093, 3.5196212224289765 … 0.9122143221397975, 0.
9122143221401331, 0.9122143221409135, 0.912214322140921, 0.9122143221408562
, 0.9122143221405419, 0.9122143221403283, 0.9122143221397983, 0.91221432213
941, 0.9122143221386628]
[2.6927687450844022, 2.692768793182424, 2.692768843173625, 2.6927688931925
88, 2.6927689412331732, 2.6927689852255514, 2.692769023140844, 2.6927690531
18851, 2.692769073610191, 2.692769083511124 … 1.1449286955907596, 1.14492
86955911178, 1.1449286955915916, 1.1449286955916638, 1.1449286955916267, 1.
1449286955914064, 1.1449286955911533, 1.144928695590697, 1.1449286955902385
, 1.1449286955895766]
[1.6638910873620831, 1.663891135459961, 1.663891185450922, 1.6638912354696
964, 1.6638912835102115, 1.6638913275024467, 1.6638913654175496, 1.66389139
53954915, 1.6638914158866707, 1.6638914257875743 … 1.6489243165014276, 1.
6489243165018714, 1.6489243165020488, 1.648924316502046, 1.6489243165021266
, 1.6489243165019023, 1.6489243165019143, 1.6489243165013963, 1.64892431650
0767, 1.6489243165002527]
[1.1678842944971668, 1.1678843425948444, 1.1678843925856235, 1.16788444260
42179, 1.1678844906445345, 1.1678845346366145, 1.1678845725515745, 1.167884
602529421, 1.1678846230205415, 1.1678846329214068 … 2.0516734038388442, 2
.0516734038391453, 2.0516734038392626, 2.0516734038392577, 2.05167340383931
45, 2.0516734038391613, 2.0516734038391777, 2.0516734038388256, 2.051673403
8383984, 2.0516734038380546]
[0.6069466972096277, 0.6069467453083509, 0.6069467952970033, 0.60694684531
44846, 0.6069468933526092, 0.6069469373461568, 0.6069469752616932, 0.606947
00523755, 0.606947025727838, 0.6069470356287272 … 2.7995778606854413, 2.7
995778606865453, 2.799577860686138, 2.799577860686269, 2.799577860686521, 2
.7995778606861785, 2.7995778606871626, 2.799577860686091, 2.799577860687871
, 2.799577860687257]
[0.47215555277672316, 0.4721556008742116, 0.47215565086411865, 0.472155700
88207404, 0.4721557489214685, 0.4721557929135511, 0.4721558308283381, 0.472
15586080552696, 0.4721558812963327, 0.47215589119712137 … 3.1676671971764
154, 3.16766719717659, 3.167667197176464, 3.1676671971764727, 3.16766719717
6529, 3.1676671971764807, 3.167667197176733, 3.1676671971765695, 3.16766719
71770274, 3.1676671971769923]
[0.4055628704418776, 0.4055629185389692, 0.405562968529316, 0.405563018547
4488, 0.40556306658728697, 0.4055631105788928, 0.4055631484934412, 0.405563
1784710878, 0.4055631989620743, 0.40556320886283254 … 3.603368151456505,
3.603368151456402, 3.6033681514563685, 3.6033681514563423, 3.60336815145633
8, 3.6033681514563765, 3.6033681514563987, 3.603368151456506, 3.60336815145
6546, 3.6033681514566753]
[0.3995929112700888, 0.3995929593673295, 0.39959300935761466, 0.3995930593
7575146, 0.39959310741552356, 0.39959315140729207, 0.39959318932193655, 0.3
995932192994949, 0.3995932397904544, 0.39959324969122934 … 3.910174923881
8027, 3.9101749238817978, 3.9101749238817582, 3.91017492388175, 3.910174923
881759, 3.910174923881765, 3.9101749238818266, 3.910174923881835, 3.9101749
238819465, 3.910174923881991]
[0.42180482514855694, 0.42180487324580856, 0.421804923236389, 0.4218049732
5473093, 0.42180502129478975, 0.4218050652865224, 0.42180510320119247, 0.42
18051331789871, 0.4218051536700632, 0.4218051635708545 … 4.32769488090108
9, 4.327694880901067, 4.327694880901113, 4.3276948809011175, 4.327694880901
103, 4.327694880901107, 4.327694880901029, 4.327694880901046, 4.32769488090
09035, 4.327694880900876]
Visualizing the solution (works best in a terminal):
@gif for t in sol1.t[1]:0.1:sol1.t[end] off = N^2 solt = sol1(t) plt1 = surface(reshape(solt[1:off], N, N), zlims=(0, 5), leg=false) surface!(plt1, reshape(solt[off+1:end], N, N), zlims=(0, 5), leg=false) display(plt1) end
Part 2: Optimizing the BRUSS Code
function brusselator_2d_loop(du, u, p, t) A, B, α, xyd = p dx = step(xyd) N = length(xyd) α = α/dx^2 limit = a -> let N=N a == N+1 ? 1 : a == 0 ? N : a end II = LinearIndices((N, N, 2)) @inbounds begin for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] i = I[1] j = I[2] ip1 = limit(i+1) im1 = limit(i-1) jp1 = limit(j+1) jm1 = limit(j-1) ii1 = II[i,j,1] ii2 = II[i,j,2] du[II[i,j,1]] = α*(u[II[im1,j,1]] + u[II[ip1,j,1]] + u[II[i,jp1,1]] + u[II[i,jm1,1]] - 4u[ii1]) + B + u[ii1]^2*u[ii2] - (A + 1)*u[ii1] + brusselator_f(x, y, t) du[II[i,j,2]] = α*(u[II[im1,j,2]] + u[II[ip1,j,2]] + u[II[i,jp1,2]] + u[II[i,jm1,2]] - 4u[II[i,j,2]]) + A*u[ii1] - u[ii1]^2*u[ii2] end end nothing end prob2 = ODEProblem(brusselator_2d_loop, u0, tspan, p) sol2 = @time solve(prob2, TRBDF2())
9.398752 seconds (12.04 M allocations: 710.365 MiB, 4.81% gc time)
sol2_2 = @time solve(prob2, CVODE_BDF())
39.956680 seconds (2.82 M allocations: 171.969 MiB, 0.05% gc time)
retcode: Success
Interpolation: 3rd order Hermite
t: 259-element Array{Float64,1}:
0.0
5.653338356947935e-11
5.65390369078363e-7
2.572950290187468e-6
4.580510211296573e-6
7.750405709122968e-6
1.4302881311757801e-5
4.0134710329210774e-5
6.596653934666375e-5
9.179836836411672e-5
⋮
20.115148615356873
20.37146777739318
20.627786939429484
20.88410610146579
21.140425263502095
21.3967444255384
21.653063587574707
21.909382749611012
22.0
u: 259-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 2.5250877095783344,
2.2620667554742258, 1.9735248771761977, 1.665005111992191, 1.34336401668221
03, 1.0172186542655526, 0.6977464117458191, 0.4003323380813969, 0.148922584
53196738, 0.0]
[6.598117452446044e-8, 6.598117452452795e-8, 6.59811745246415e-8, 6.598117
452477602e-8, 6.59811745249205e-8, 6.598117452506815e-8, 6.598117452521404e
-8, 6.598117452535298e-8, 6.598117452548406e-8, 6.598117452560301e-8 … 2.
52508769326683, 2.26206674160907, 1.9735248663224956, 1.6650051048635588, 1
.3433640142351067, 1.0172186578909652, 0.6977464237297119, 0.40033236307490
11, 0.1489226402117897, 8.090751419307918e-8]
[0.0006587624140132232, 0.0006587624209141625, 0.0006587624326779875, 0.00
06587624465990196, 0.000658762461542147, 0.000658762476747418, 0.0006587624
918359335, 0.000658762506258849, 0.0006587625198234356, 0.00065876253213972
14 … 2.524924606099231, 2.261928123190901, 1.9734163705350949, 1.66493387
301444, 1.343339621372235, 1.0172550391384978, 0.697866520902045, 0.4005832
5044085197, 0.14947914535482773, 0.0008077938023952315]
[0.0029885797701746284, 0.002988579862875722, 0.002988580024412284, 0.0029
885802153883506, 0.0029885804198764242, 0.002988580626391515, 0.00298858083
42275124, 0.0029885810313851017, 0.0029885812191603584, 0.00298858138811432
6 … 2.5243457007698074, 2.261436153310274, 1.9730314375969087, 1.66468138
01159555, 1.3432536595673659, 1.0173853107174544, 0.6982952263862713, 0.401
48198029752397, 0.1514529070863203, 0.003664729223250155]
[0.005306141508030984, 0.005306142690205768, 0.005306142257972833, 0.00530
6142797468204, 0.005306143401977575, 0.00530614493691719, 0.005306148087986
683, 0.0053061452495214295, 0.005306149230355061, 0.005306149723653397 …
2.5237671078891086, 2.2609445550730123, 1.972646967738447, 1.66442950038285
24, 1.3431685811161493, 1.0175170325180882, 0.698726972801953, 0.4023908498
594301, 0.15342343126294944, 0.0065066953025461555]
[0.008938531142281401, 0.00893853412998814, 0.008938533324753955, 0.008938
534851373863, 0.008938536624027108, 0.008938540672423655, 0.008938548169476
477, 0.00893854209409423, 0.008938551523646868, 0.008938552952542707 … 2.
5228542226523203, 2.2601691664198333, 1.9720409392385245, 1.664033168895644
7, 1.3430362416837864, 1.0177283298803772, 0.6994157785454703, 0.4038476392
5678996, 0.15652765095572066, 0.010961135026978663]
[0.016344145375892745, 0.016344150679504905, 0.016344153382255998, 0.01634
4159160398462, 0.01634416578191356, 0.016344175003834897, 0.016344186710899
49, 0.016344185944460035, 0.01634419915335025, 0.01634420450260158 … 2.52
09699686322784, 2.2585696653032357, 1.9707923492854083, 1.6632194124484905,
1.3427706997876703, 1.0181786440181924, 0.7008691316897557, 0.406937430466
94203, 0.16291636532495113, 0.020043247383951576]
[0.04433812091521405, 0.04433815867165508, 0.04433820623536491, 0.04433827
157330291, 0.044338347135081634, 0.044338431939997466, 0.044338515477739565
, 0.04433857623295058, 0.044338655086562505, 0.04433871696545116 … 2.5135
78240692516, 2.252307885072983, 1.9659255801385158, 1.66008673108698, 1.341
8375910046947, 1.0201544536483305, 0.7070145269795118, 0.4198921130658216,
0.18773948635125548, 0.054381120820077296]
[0.07058737985116977, 0.07058664484029947, 0.0705875413280103, 0.070587855
81917236, 0.0705879729541195, 0.070587289183245, 0.07058596349206522, 0.070
58829612748044, 0.07058625771406031, 0.0705864535189533 … 2.5062449271565
765, 2.2461167960972164, 1.961149034074087, 1.6570783972413137, 1.341096371
9495257, 1.0224707341431893, 0.7138163912529947, 0.4338374639900661, 0.2119
8260672435237, 0.08658995681881776]
[0.09534502498283606, 0.09534358538038319, 0.09534532706763858, 0.09534608
102785042, 0.09534625281411645, 0.09534489421692781, 0.09534277774807996, 0
.0953468096914476, 0.09534337195439287, 0.09534381035303054 … 2.498970563
5047277, 2.2399977238708053, 1.9564660515271148, 1.6542026515183286, 1.3405
654586114462, 1.0251552440356688, 0.7212413308882024, 0.44843068771835926,
0.23564736424496766, 0.1169787925933724]
⋮
[0.6323810451914815, 0.632427052649651, 0.6324748742136165, 0.632522725663
7151, 0.632568687588502, 0.6326107791949938, 0.6326470583214443, 0.63267574
42982498, 0.6326953532738601, 0.632704828378052 … 2.770008193556327, 2.77
0008208734769, 2.7700082189120816, 2.7700082237279133, 2.7700082230441634,
2.7700082168913203, 2.7700082054620894, 2.770008189124797, 2.77000816842502
84, 2.7700081440708226]
[0.5226275087207661, 0.5226734342910327, 0.5227211719697351, 0.52276894073
76451, 0.5228148244346382, 0.5228568454418254, 0.5228930645349704, 0.522921
7036079294, 0.5229412808642269, 0.5229507408309406 … 3.0270587055167733,
3.0270586753006343, 3.0270586550848355, 3.027058645472268, 3.02705864680864
9, 3.027058659063031, 3.0270586818168277, 3.0270587143052663, 3.02705875544
77606, 3.0270588038594153]
[0.4579604087410305, 0.458006296802263, 0.4580539960946591, 0.458101726935
9897, 0.45814757462705746, 0.458189563162972, 0.45822575460980086, 0.458254
3720215345, 0.4582739347433197, 0.458283387662069 … 3.2613243481002447, 3
.2613242970972784, 3.261324263504071, 3.2613242471006543, 3.261324249141863
7, 3.261324269879769, 3.261324308306207, 3.2613243629016937, 3.261324432013
5792, 3.2613245136746594]
[0.4243249949117141, 0.4243708716334364, 0.4244185592849013, 0.42446627872
460485, 0.42451211572444175, 0.42455409459691923, 0.42459027785523756, 0.42
46188889172981, 0.42463844725967403, 0.42464789817336857 … 3.476979389713
562, 3.4769793324322897, 3.476979294466571, 3.476979276125042, 3.4769792785
14862, 3.476979301782493, 3.4769793449316233, 3.476979406359194, 3.47697948
4132899, 3.4769795758778326]
[0.4100847046913288, 0.4101305921364255, 0.4101782907401593, 0.41022602107
792533, 0.41027186852471176, 0.4103138568061762, 0.41035004810598086, 0.410
37866551793173, 0.4103982280350711, 0.4104076810589709 … 3.67907975182885
67, 3.6790797004957008, 3.679079666185172, 3.679079649844482, 3.67907965210
2726, 3.679079672926244, 3.679079711585095, 3.6790797667656774, 3.679079836
642917, 3.679079918887027]
[0.4073307682105559, 0.4073766856477378, 0.4074244149682815, 0.40747217562
354593, 0.4075180518070638, 0.40756006601996625, 0.40759627938780474, 0.407
62491406752893, 0.4076444882216691, 0.407653946842301 … 3.87153703747294,
3.8715370027640295, 3.8715369794090573, 3.871536968411483, 3.8715369700022
073, 3.8715369840696034, 3.8715370102065285, 3.871537047588856, 3.871537094
929266, 3.871537150541017]
[0.4116592888096325, 0.41170525037260386, 0.4117530249037378, 0.4118008301
185253, 0.41184674845518116, 0.41188880071685885, 0.4119250464392734, 0.411
9537063944312, 0.41197329765271995, 0.4119827644270476 … 4.05641774845540
05, 4.05641773820175, 4.056417731164946, 4.056417727958729, 4.0564177284843
43, 4.056417732631651, 4.056417740353393, 4.056417751458394, 4.056417765516
157, 4.056417781927513]
[0.421178883717942, 0.42122490081186037, 0.42127273224506245, 0.4213205935
2132606, 0.42136656486420737, 0.42140866497435914, 0.421444951375888, 0.421
4736430958104, 0.42149325587312986, 0.4215027328817772 … 4.23412681426794
4, 4.234126834790803, 4.234126848384466, 4.234126854951361, 4.2341268540976
16, 4.234126845769279, 4.234126830318249, 4.234126808309182, 4.234126780420
57, 4.234126747488131]
[0.4255569709603833, 0.42560301293838904, 0.42565086989133055, 0.425698756
25778114, 0.4257447512526154, 0.4257868727264482, 0.42582317727042507, 0.42
585188312176897, 0.42587150553824715, 0.42588098706735406 … 4.29522520924
0007, 4.295225243549162, 4.295225266516079, 4.295225277416021, 4.2952252758
90809, 4.295225261987599, 4.295225236161456, 4.2952251992579855, 4.29522515
2494508, 4.295225097441233]
Part 3: Exploiting Jacobian Sparsity with Color Differentiation
using SparseDiffTools, SparsityDetection sparsity_pattern = jacobian_sparsity(brusselator_2d_loop,similar(u0),u0,p,2.0)
Explored path: SparsityDetection.Path(Bool[], 1)
jac_sp = sparse(sparsity_pattern) jac = Float64.(jac_sp) colors = matrix_colors(jac) prob3 = ODEProblem(ODEFunction(brusselator_2d_loop, colorvec=colors,jac_prototype=jac_sp), u0, tspan, p) sol3 = @time solve(prob3, TRBDF2())
ERROR: InexactError: Bool(-38444.40000000001)
(Optional) Part 4: Structured Jacobians
(Optional) Part 5: Automatic Symbolicification and Analytical Jacobian
Part 6: Utilizing Preconditioned-GMRES Linear Solvers
using DiffEqOperators using Sundials using AlgebraicMultigrid: ruge_stuben, aspreconditioner, smoothed_aggregation prob6 = ODEProblem(ODEFunction(brusselator_2d_loop, jac_prototype=JacVecOperator{Float64}(brusselator_2d_loop, u0)), u0, tspan, p) II = Matrix{Float64}(I, N, N) Op = kron(Matrix{Float64}(I, 2, 2), kron(D2, II) + kron(II, D2)) Wapprox = -I+Op #ml = ruge_stuben(Wapprox) ml = smoothed_aggregation(Wapprox) precond = aspreconditioner(ml) sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES())); # no preconditioner
12.849774 seconds (17.64 M allocations: 5.066 GiB, 2.08% gc time)
sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES(Pl=lu(Wapprox)))); # sparse LU
3.459602 seconds (6.69 M allocations: 680.707 MiB, 1.93% gc time)
sol_trbdf2 = @time solve(prob6, TRBDF2(linsolve=LinSolveGMRES(Pl=precond))); # AMG
3.452518 seconds (7.46 M allocations: 535.606 MiB, 1.82% gc time)
sol_cvodebdf = @time solve(prob2, CVODE_BDF(linear_solver=:GMRES));
1.641359 seconds (2.05 M allocations: 203.154 MiB, 1.15% gc time)
retcode: Success
Interpolation: 3rd order Hermite
t: 1928-element Array{Float64,1}:
0.0
5.653338356947935e-11
5.65390369078363e-7
2.57289768832925e-6
4.5804050075801365e-6
7.750237020809387e-6
1.430297275427241e-5
3.947737124517274e-5
6.465176973607307e-5
8.982616822697341e-5
⋮
21.74842347830861
21.75377749261479
21.770859233529922
21.787940974445053
21.800484442937336
21.813027911429618
21.844693963335253
21.950342098966686
22.0
u: 1928-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 2.5250877095783344,
2.2620667554742258, 1.9735248771761977, 1.665005111992191, 1.34336401668221
03, 1.0172186542655526, 0.6977464117458191, 0.4003323380813969, 0.148922584
53196738, 0.0]
[6.598117452445592e-8, 6.598117452452327e-8, 6.598117452463697e-8, 6.59811
7452477148e-8, 6.598117452491596e-8, 6.598117452506346e-8, 6.59811745252089
2e-8, 6.598117452534844e-8, 6.598117452547893e-8, 6.598117452559789e-8 …
2.52508769326683, 2.26206674160907, 1.9735248663224956, 1.6650051048635588,
1.3433640142351067, 1.0172186578909652, 0.6977464237297119, 0.400332363074
9011, 0.1489226402117897, 8.090751419307378e-8]
[0.0006587624172964451, 0.0006587624240399749, 0.0006587624354243716, 0.00
06587624488920072, 0.0006587624633585051, 0.0006587624781271849, 0.00065876
24926919035, 0.000658762506662457, 0.0006587625197283628, 0.000658762531638
6175 … 2.5249246059817017, 2.2619281230446204, 1.973416370341658, 1.66493
38727357148, 1.3433396209127182, 1.0172550381632066, 0.69786651676992, 0.40
05832569059138, 0.14947914564976048, 0.0008077937954087897]
[0.0029885190994651673, 0.0029885191897600113, 0.0029885193423737177, 0.00
29885195232632685, 0.002988519717653178, 0.002988519916140473, 0.0029885201
11906188, 0.002988520299697163, 0.0029885204753350315, 0.002988520635442596
3 … 2.524345713343465, 2.261436162970071, 1.973031443407278, 1.6646813805
594152, 1.343253651586457, 1.017385285060445, 0.698295128369312, 0.40148209
36914462, 0.15145286366257169, 0.0036646545802277134]
[0.005306022356073513, 0.005306022621712586, 0.0053060230717453995, 0.0053
06023607201543, 0.005306024183137538, 0.005306024771435465, 0.0053060253517
8454, 0.005306025908561257, 0.005306026429347995, 0.005306026904112427 …
2.5237671339498102, 2.260944575638989, 1.9726469809636227, 1.66442950338789
62, 1.343168568455559, 1.0175169867522391, 0.6987268229191962, 0.4023910083
0343055, 0.15342334401815053, 0.0065065496316772335]
[0.00893834350995413, 0.00893834427987621, 0.00893834558850756, 0.00893834
7153465821, 0.008938348839052852, 0.008938350561768801, 0.00893835226170432
, 0.008938353892881468, 0.008938355418797277, 0.008938356809979058 … 2.52
28542649298107, 2.2601692000771423, 1.9720409612790466, 1.66403317465477, 1
.3430362221641639, 1.017728256881341, 0.6994155885173241, 0.403847823621840
9, 0.15652751916175822, 0.01096090625912015]
[0.01634429202116057, 0.01634429489952914, 0.016344299822335723, 0.0163443
0576423857, 0.016344312182172048, 0.016344318748721217, 0.01634432523215556
5, 0.01634433145550153, 0.01634433727856935, 0.016344342588309163 … 2.520
9699337854867, 2.258569632374003, 1.970792317068139, 1.6632193768804995, 1.
3427706450978134, 1.0181785268879573, 0.7008693575988351, 0.406937285311047
1, 0.16291664462188093, 0.0200434281564615]
[0.043646263363862924, 0.043646292669642976, 0.04364634410002619, 0.043646
408299418085, 0.043646478690166544, 0.0436465511401205, 0.04364662287963492
, 0.04364669185715508, 0.043646756468291935, 0.04364681542761546 … 2.5137
65529699384, 2.2524662162038327, 1.96604803719552, 1.6601642835144645, 1.34
18578790218518, 1.0201021423806298, 0.7068529962881926, 0.41955020644048985
, 0.1871136496191611, 0.05353236358962729]
[0.0692823968869293, 0.06928248583818256, 0.06928264396753456, 0.069282844
55502629, 0.06928306665139533, 0.06928329638678168, 0.06928352441574115, 0.
06928374395054913, 0.06928394975415433, 0.06928413765698549 … 2.506616550
182622, 2.2464297079132223, 1.9613890075261213, 1.6572267432709649, 1.34112
728871286, 1.0223507377796497, 0.7134630924575552, 0.43311194946272524, 0.2
1076011686633297, 0.08498786962833081]
[0.09349042837387436, 0.09349061562324566, 0.09349095207267255, 0.09349138
44519549, 0.09349186761490513, 0.0934923699926786, 0.09349286998656095, 0.0
9349335206219281, 0.09349380438256111, 0.09349421759820462 … 2.4995234731
16373, 2.240461415578458, 1.956818523753225, 1.65441487227758, 1.3405970544
326884, 1.0249482572805424, 0.720666723995856, 0.447302887948819, 0.2338561
2858767357, 0.11470081104803237]
⋮
[0.3797955575005551, 0.3798129120680182, 0.37983160912795294, 0.3798515167
9769735, 0.3798700744173499, 0.3798877958739413, 0.37990413586835, 0.379917
73881418767, 0.3799267950426437, 0.3799309569353526 … 3.4574052823485664,
3.4573936286619094, 3.457386402854468, 3.4573827565631308, 3.4573831040940
81, 3.4573876418782974, 3.457395923872096, 3.4574086494131495, 3.4574240121
578605, 3.4574420811420903]
[0.37958401654816143, 0.3796201588394765, 0.37965786768047793, 0.379695834
60132494, 0.3797327833587292, 0.379767041189448, 0.37979645035002474, 0.379
81940648060253, 0.37983513643422023, 0.3798427397648622 … 3.4617077059990
3, 3.4616995134797413, 3.4616938254044576, 3.4616910792093196, 3.4616914818
578866, 3.4616948057305494, 3.46170100881632, 3.4617097544879343, 3.4617206
586988303, 3.4617334162991207]
[0.37898131515329514, 0.37904455623143324, 0.3791102395621268, 0.379171848
92573214, 0.3792278014621423, 0.37927854226433344, 0.3793222968917526, 0.37
935935408014265, 0.3793873038928219, 0.37940252613199893 … 3.475400682841
286, 3.4754021781703925, 3.475405250938518, 3.4754068566952996, 3.475406283
915746, 3.475404220521274, 3.475401075649539, 3.4754010132857456, 3.4754017
445012275, 3.4754014416136703]
[0.3784028728681648, 0.37847935063114085, 0.37855821248433963, 0.378635433
8343628, 0.3787103744504435, 0.37877900390104313, 0.3788373625105384, 0.378
88198986151156, 0.37891190751041437, 0.37892669234514176 … 3.489080362945
1366, 3.489089294219412, 3.489094650501359, 3.489097497162135, 3.4890973438
878747, 3.489093753047733, 3.489087808019096, 3.4890784385409153, 3.4890666
40825046, 3.489052679244454]
[0.37802806562457314, 0.3781081068261598, 0.37819023494543214, 0.378268702
7300955, 0.37834355080180027, 0.37841307565234206, 0.37847181048385786, 0.3
785172210264568, 0.37854884544138667, 0.3785659311574103 … 3.499071751303
758, 3.4990834581052974, 3.4990920535881895, 3.499096339156102, 3.499095815
52635, 3.4990905269259645, 3.49908103431579, 3.4990696874410214, 3.49905517
5741332, 3.4990372237525125]
[0.3776909382220414, 0.37776351115109563, 0.37783859475543263, 0.377912506
00268866, 0.37798520447195416, 0.37805140767486484, 0.37810798564087067, 0.
37815118643386414, 0.3781792196305945, 0.37819286030587723 … 3.5090277181
657115, 3.5090394826004574, 3.5090464815977596, 3.509050187200965, 3.509049
8417172786, 3.509045416275706, 3.5090377144093337, 3.509024825081179, 3.509
0089285477553, 3.508989993790065]
[0.3768537939605036, 0.376915458761548, 0.3769793531358606, 0.377042729545
79915, 0.37710281438118, 0.3771570791988641, 0.37720388440231123, 0.3772415
338113573, 0.37726733457251116, 0.3772798488606232 … 3.5341352486353617,
3.5341466559569223, 3.534154680705109, 3.534158571433794, 3.534157981996896
2, 3.5341533014526756, 3.5341445069751627, 3.534132413334222, 3.53411749417
5883, 3.5340999041834755]
[0.3740747443719409, 0.3741028206103584, 0.3741325900228058, 0.37416663902
824293, 0.3742006984167679, 0.3742310255447786, 0.3742579572053362, 0.37427
89270608717, 0.3742919258904212, 0.37429598291067795 … 3.6178769692717765
, 3.6178851739555844, 3.617889139403817, 3.6178908996243493, 3.617890994937
32, 3.617888791395569, 3.6178843129218543, 3.617874139997575, 3.61786105897
74163, 3.617847139077364]
[0.37278515952965974, 0.3728179321140629, 0.3728522412070423, 0.3728876482
7093895, 0.37291708637372784, 0.37294295966405455, 0.37296625674297096, 0.3
7298877738348996, 0.3730068685691214, 0.3730153360421533 … 3.657191327231
621, 3.6571972377609487, 3.65720344716516, 3.6572059278358457, 3.6572051957
615703, 3.6572021410971995, 3.6571953691210877, 3.657189705548747, 3.657182
8565260756, 3.6571749711705777]
Part 7: Exploring IMEX and Exponential Integrator Techniques (E)
function laplacian2d(du, u, p, t) A, B, α, xyd = p dx = step(xyd) N = length(xyd) du = reshape(du, N, N, 2) u = reshape(u, N, N, 2) @inbounds begin α = α/dx^2 limit = a -> let N=N a == N+1 ? 1 : a == 0 ? N : a end for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] i = I[1] j = I[2] ip1 = limit(i+1) im1 = limit(i-1) jp1 = limit(j+1) jm1 = limit(j-1) du[i,j,1] = α*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) du[i,j,2] = α*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) end end nothing end function brusselator_reaction(du, u, p, t) A, B, α, xyd = p dx = step(xyd) N = length(xyd) du = reshape(du, N, N, 2) u = reshape(u, N, N, 2) @inbounds begin for I in CartesianIndices((N, N)) x = xyd[I[1]] y = xyd[I[2]] i = I[1] j = I[2] du[i,j,1] = B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t) du[i,j,2] = A*u[i,j,1] - u[i,j,1]^2*u[i,j,2] end end nothing end prob7 = SplitODEProblem(laplacian2d, brusselator_reaction, u0, tspan, p) sol7 = @time solve(prob7, KenCarp4())
7.972783 seconds (24.90 M allocations: 1.342 GiB, 1.55% gc time)
M = MatrixFreeOperator((du,u,p)->laplacian2d(du, u, p, 0), (p,), size=(2*N^2, 2*N^2), opnorm=1000) prob7_2 = SplitODEProblem(M, brusselator_reaction, u0, tspan, p) sol7_2 = @time solve(prob7_2, ETDRK4(krylov=true), dt=1)
17.154728 seconds (63.51 M allocations: 3.244 GiB, 4.70% gc time)
prob7_3 = SplitODEProblem(DiffEqArrayOperator(Op), brusselator_reaction, u0, tspan, p) sol7_3 = solve(prob7_3, KenCarp4());
retcode: Success
Interpolation: 3rd order Hermite
t: 411-element Array{Float64,1}:
0.0
2.888065882181688e-5
0.00031768724703998567
0.0027059840153272462
0.008731809471821457
0.02179419538626385
0.046044933347512154
0.08666612492700337
0.14668700100553306
0.2429330434920095
⋮
21.642678211001897
21.669167858631038
21.69565750626018
21.745650901936376
21.842907495920542
21.862544058930474
21.883418684049154
21.924650607598824
22.0
u: 411-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 2.5250877095783344,
2.2620667554742258, 1.9735248771761977, 1.665005111992191, 1.34336401668221
03, 1.0172186542655526, 0.6977464117458191, 0.4003323380813969, 0.148922584
53196738, 0.0]
[0.00338181359279749, 0.003381814525055848, 0.0033818161026810765, 0.00338
18179763281625, 0.0033818199908016255, 0.003381822048157553, 0.003381824077
5272284, 0.0033818260243597817, 0.00338182784528272, 0.003381829505247768
… 2.524255066788988, 2.261359157113013, 1.9729712293922774, 1.664641937979
7444, 1.3432403152050174, 1.0174058169618412, 0.6983623120712237, 0.4016217
8277595556, 0.1517635265052733, 0.004115818959992079]
[0.03568250376202336, 0.03568266931785142, 0.035682955832227074, 0.0356833
0729754052, 0.03568368941345571, 0.03568408140938937, 0.03568446895364048,
0.035684841242777064, 0.035685189766591796, 0.03568550768135772 … 2.51599
73402014577, 2.2543598849568935, 1.967523764161279, 1.6611194391842758, 1.3
421496263916268, 1.0195052150551678, 0.704979868148862, 0.4156649097840782,
0.17979934864987152, 0.04350153844352491]
[0.2388870086712341, 0.23890840746265182, 0.23894852951279852, 0.239002926
71217413, 0.23906668747399049, 0.23913531679710034, 0.23920518662122434, 0.
2392735177564327, 0.2393382102869443, 0.23939766131788623 … 2.45177455071
75612, 2.2013122441306154, 1.9285687665074842, 1.6402584486989873, 1.345213
7788742935, 1.0552180035855303, 0.7857635817095209, 0.5558670473603251, 0.3
8607263451905216, 0.29539374940409546]
[0.5564106150796982, 0.5566108464330898, 0.5569986710045511, 0.55755022257
69137, 0.5582336728011149, 0.5590119152749861, 0.559846526329856, 0.5606997
95457231, 0.561537663173318, 0.5623298108904917 … 2.3212076996976343, 2.1
072090690298686, 1.880355997024359, 1.6487387710243555, 1.4216785859442722,
1.2094079500272332, 1.0225411460886549, 0.8713229853793858, 0.764708535522
895, 0.7095402391463423]
[0.9342714734084758, 0.9349598918146564, 0.936310382519549, 0.938262165081
3462, 0.9407415908113002, 0.9436441970638244, 0.9468551430516244, 0.9502481
957941415, 0.9536902430916663, 0.957050094472094 … 2.1614414096004184, 2.
027111451991201, 1.8887980771792972, 1.751800690784729, 1.6215018392791603,
1.5031551617020042, 1.4016704666279294, 1.321363816947498, 1.2657330373145
712, 1.2372638912267229]
[1.2566772231894945, 1.2576431810047317, 1.2595383902093107, 1.26229311347
9925, 1.2658028435853461, 1.2699372640323565, 1.274540022703064, 1.27943691
0349005, 1.2844410563581503, 1.2893600761202952 … 2.0368655159549163, 1.9
851294449856716, 1.9322949399059546, 1.8803923426225904, 1.8314233258837715
, 1.7872824916440149, 1.749687377223869, 1.7201110332114784, 1.699714353288
5108, 1.689307976958115]
[1.4173304157928486, 1.4178104626035777, 1.4187536323652747, 1.42012345899
5461, 1.421869253600005, 1.4239236083651097, 1.4262077943444462, 1.42863463
98967849, 1.4311109979822425, 1.4335434070534994 … 1.984707244631404, 1.9
744843346627017, 1.9640642411192235, 1.9538390653254547, 1.9442016835060763
, 1.9355229382392862, 1.9281377603979366, 1.9223312732585927, 1.91833313347
7084, 1.9162957064303392]
[1.4079292375651262, 1.4080181845337447, 1.4081926434529766, 1.40844579737
3808, 1.40876792341743, 1.4091467072462807, 1.4095677840064875, 1.410015143
058539, 1.410471755723265, 1.4109202040012978 … 2.006848717030137, 2.0058
633341206855, 2.004859241773954, 2.003874429852624, 2.0029462538824254, 2.0
02109807224481, 2.0013971225351015, 2.000835759821026, 2.000448425718132, 2
.0002507800301195]
[1.3008722049712922, 1.3008801828116636, 1.3008945778368928, 1.30091484490
6319, 1.3009396166641172, 1.3009679104321659, 1.3009992327229458, 1.3010331
324351492, 1.3010682282392434, 1.3011028322919374 … 2.084972574467494, 2.
0849211490150776, 2.0848702925700606, 2.0848214163333547, 2.08477743413165,
2.084738198481765, 2.0847036216234702, 2.084673859207271, 2.08465246394274
, 2.0846409823084646]
⋮
[0.39890859367367754, 0.3993441673525077, 0.3997975806788334, 0.4002518902
321085, 0.4006888182480096, 0.4010893923770546, 0.40143496888258334, 0.4017
0842611233254, 0.401895469313635, 0.40198591841037423 … 3.974042467040040
7, 3.9740400700975256, 3.974038463649698, 3.974037701900219, 3.974037809267
5783, 3.9740387819095946, 3.974040587360848, 3.974043164462041, 3.974046422
7965795, 3.974050246993257]
[0.399359295339172, 0.39974889839511, 0.40014696627279867, 0.4005384174406
543, 0.4009068845049872, 0.40123752344283037, 0.40151687180652956, 0.401733
60655661844, 0.401879356549643, 0.4019484860747507 … 3.9932059688091033,
3.9932038478068232, 3.9932024382289315, 3.9932017744025847, 3.9932018738678
297, 3.99320272691542, 3.993204321872256, 3.993206602084878, 3.993209523141
2905, 3.9932129859626735]
[0.3993686665009743, 0.39980665253816694, 0.4002705836448292, 0.4007431063
3130167, 0.40120406441906215, 0.4016317190585153, 0.4020041495019104, 0.402
30097607591025, 0.40250515220793087, 0.4026046039165201 … 4.0122858125097
78, 4.012284245475603, 4.012283191763723, 4.01228269079853, 4.0122827626693
17, 4.012283404160809, 4.012284589905506, 4.012286273787084, 4.012288391829
275, 4.012290860181647]
[0.4010117037388889, 0.4014585319070676, 0.40192471908753297, 0.4023929107
571625, 0.40284400510226775, 0.4032581888842881, 0.4036158867371329, 0.4038
991410242888, 0.4040929782391569, 0.40418678544141134 … 4.048074184182981
5, 4.048072793476128, 4.048071860809855, 4.0480714180078925, 4.048071480580
622, 4.04807204591452, 4.048073094164494, 4.0480745886306035, 4.04807647493
2387, 4.048078682581573]
[0.4040391360250253, 0.4044781455356283, 0.40493508964458735, 0.4053928897
130187, 0.40583313109968033, 0.4062367046948323, 0.4065848385324651, 0.4068
602985846021, 0.40704869896301227, 0.40713979724380694 … 4.11688589803744
4, 4.116885469878567, 4.116885180300814, 4.116885042860157, 4.1168850632805
38, 4.1168852404840575, 4.116885566552351, 4.116886026211559, 4.11688659676
3158, 4.116887250002496]
[0.4047700956071005, 0.40519815829045214, 0.4056399542593235, 0.4060787935
4982956, 0.40649646573647547, 0.4068752731237485, 0.4071985546331393, 0.407
45169655229097, 0.40762327673507454, 0.4077054230583039 … 4.1306436987736
2, 4.130643465603472, 4.130643313409505, 4.130643244570406, 4.1306432586828
51, 4.130643353557858, 4.130643529891911, 4.130643781395985, 4.130644101166
296, 4.130644470170689]
[0.405168629966455, 0.4055798821600529, 0.40601390501288664, 0.40645506728
755565, 0.4068852428751017, 0.4072847292800238, 0.4076332639179229, 0.40791
16747690355, 0.40810359543776764, 0.40819724273804125 … 4.145219717875681
, 4.145219791546891, 4.145219834102426, 4.14521985299266, 4.145219852795447
, 4.145219832891524, 4.145219787681931, 4.145219708636722, 4.14521958516251
4, 4.145219403944229]
[0.40704121426682677, 0.40748978803601077, 0.4079566870477457, 0.408424669
7332447, 0.40887498057174415, 0.4092881700600071, 0.40964487371286185, 0.40
9927358248609, 0.4101206777465779, 0.41021423546960456 … 4.17385578570464
1, 4.17385633196801, 4.173856695023444, 4.173856867186141, 4.17385684409214
35, 4.173856625935988, 4.173856217791514, 4.173855628428838, 4.173854869192
83, 4.173853955130709]
[0.4104280757891703, 0.4108700639501849, 0.41132984460416, 0.4117903516752
1066, 0.41223306465343723, 0.41263893115959094, 0.41298899290397084, 0.4132
6603039541165, 0.4134554904299645, 0.4135471107118313 … 4.225664869407304
, 4.225666277355228, 4.225667215715375, 4.22566766028349, 4.225667599839396
, 4.225667035917906, 4.225665982952057, 4.225664468355488, 4.22566253340619
55, 4.225660233615377]
Part 8: Work-Precision Diagrams for Benchmarking Solver Choices
using DiffEqDevTools abstols = 0.1 .^ (5:8) reltols = 0.1 .^ (1:4) sol = solve(prob3,CVODE_BDF(linear_solver=:GMRES),abstol=1/10^7,reltol=1/10^10) test_sol = TestSolution(sol) probs = [prob2, prob3, prob6] setups = [Dict(:alg=>CVODE_BDF(),:prob_choice => 1), Dict(:alg=>CVODE_BDF(linear_solver=:GMRES), :prob_choice => 1), Dict(:alg=>TRBDF2(), :prob_choice => 1), Dict(:alg=>TRBDF2(linsolve=LinSolveGMRES(Pl=precond)), :prob_choice => 3), Dict(:alg=>TRBDF2(), :prob_choice => 2) ] labels = ["CVODE_BDF (dense)" "CVODE_BDF (GMRES)" "TRBDF2 (dense)" "TRBDF2 (sparse)" "TRBDF2 (GMRES)"] wp = WorkPrecisionSet(probs,abstols,reltols,setups;appxsol=[test_sol,test_sol,test_sol],save_everystep=false,numruns=3, names=labels, print_names=true, seconds=0.5)
CVODE_BDF (dense) CVODE_BDF (GMRES) TRBDF2 (dense) TRBDF2 (sparse) TRBDF2 (GMRES)
ERROR: InexactError: Bool(-38444.40000000001)
plot(wp)
ERROR: UndefVarError: wp not defined
Part 9: GPU-Parallelism for PDEs (E)
Part 10: Adjoint Sensitivity Analysis for Gradients of PDEs
Problem 5: Global Parameter Sensitivity and Optimality with GPU and Distributed Ensembles (B)
Part 1: Implementing the Henon-Heiles System (B)
function henon(dz,z,p,t) p₁, p₂, q₁, q₂ = z[1], z[2], z[3], z[4] dp₁ = -q₁*(1 + 2q₂) dp₂ = -q₂-(q₁^2 - q₂^2) dq₁ = p₁ dq₂ = p₂ dz .= [dp₁, dp₂, dq₁, dq₂] return nothing end u₀ = [0.1, 0.0, 0.0, 0.5] prob = ODEProblem(henon, u₀, (0., 1000.)) sol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14) plot(sol, vars=[(3,4,1)], tspan=(0,100))
(Optional) Part 2: Alternative Dynamical Implmentations of Henon-Heiles (B)
function henon(ddz,dz,z,p,t) p₁, p₂ = dz[1], dz[2] q₁, q₂ = z[1], z[2] ddq₁ = -q₁*(1 + 2q₂) ddq₂ = -q₂-(q₁^2 - q₂^2) ddz .= [ddq₁, ddq₂] end p₀ = u₀[1:2] q₀ = u₀[3:4] prob2 = SecondOrderODEProblem(henon, p₀, q₀, (0., 1000.)) sol = solve(prob2, DPRKN6(), abstol=1e-10, reltol=1e-10) plot(sol, vars=[(3,4)], tspan=(0,100)) H(p, q, params) = 1/2 * (p[1]^2 + p[2]^2) + 1/2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2] - 2/3*q[2]^3) prob3 = HamiltonianProblem(H, p₀, q₀, (0., 1000.)) sol = solve(prob3, DPRKN6(), abstol=1e-10, reltol=1e-10) plot(sol, vars=[(3,4)], tspan=(0,100))
Part 3: Parallelized Ensemble Solving
In order to solve with an ensamble we need some initial conditions.
function generate_ics(E,n) # The hardcoded values bellow can be estimated by looking at the # figures in the Henon-Heiles 1964 article qrange = range(-0.4, stop = 1.0, length = n) prange = range(-0.5, stop = 0.5, length = n) z0 = Vector{Vector{typeof(E)}}() for q in qrange V = H([0,0],[0,q],nothing) V ≥ E && continue for p in prange T = 1/2*p^2 T + V ≥ E && continue z = [√(2(E-V-T)), p, 0, q] push!(z0, z) end end return z0 end z0 = generate_ics(0.125, 10) function prob_func(prob,i,repeat) @. prob.u0 = z0[i] prob end ensprob = EnsembleProblem(prob, prob_func=prob_func) sim = solve(ensprob, Vern9(), EnsembleThreads(), trajectories=length(z0)) plot(sim, vars=(3,4), tspan=(0,10))
Part 4: Parallelized GPU Ensemble Solving
In order to use GPU parallelization we must make all inputs (initial conditions, tspan, etc.) Float32 and the function definition should be in the in-place form, avoid bound checking and return nothing.
using DiffEqGPU function henon_gpu(dz,z,p,t) @inbounds begin dz[1] = -z[3]*(1 + 2z[4]) dz[2] = -z[4]-(z[3]^2 - z[4]^2) dz[3] = z[1] dz[4] = z[2] end return nothing end z0 = generate_ics(0.125f0, 50) prob_gpu = ODEProblem(henon_gpu, Float32.(u₀), (0.f0, 1000.f0)) ensprob = EnsembleProblem(prob_gpu, prob_func=prob_func) sim = solve(ensprob, Tsit5(), EnsembleGPUArray(), trajectories=length(z0))
EnsembleSolution Solution of length 1440 with uType:
DiffEqBase.ODESolution{Float32,2,Array{SubArray{Float32,1,Array{Float32,2},
Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1},Nothing,Nothing,Array{F
loat32,1},Nothing,DiffEqBase.ODEProblem{Array{Float32,1},Tuple{Float32,Floa
t32},true,DiffEqBase.NullParameters,DiffEqBase.ODEFunction{true,typeof(Main
.##WeaveSandBox#253.henon_gpu),LinearAlgebra.UniformScaling{Bool},Nothing,N
othing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Noth
ing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tup
le{}}},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsit5,DiffEqBase.Linea
rInterpolation{Array{Float32,1},Array{SubArray{Float32,1,Array{Float32,2},T
uple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}},DiffEqBase.DEStats}
Problem 6: Training Neural Stochastic Differential Equations with GPU acceleration (I)
Part 1: Constructing and Training a Basic Neural ODE
Part 2: GPU-accelerating the Neural ODE Process
Part 3: Defining and Training a Mixed Neural ODE
Part 4: Constructing a Basic Neural SDE
Part 5: Optimizing the training behavior with minibatching (E)
Information on the Build
Appendix
This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl
To locally run this tutorial, do the following commands:
using DiffEqTutorials
DiffEqTutorials.weave_file("exercises","02-workshop_solutions.jmd")Computer Information:
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqTutorials.jl/.julia
JULIA_CUDA_MEMORY_LIMIT = 536870912
JULIA_PROJECT = @.
JULIA_NUM_THREADS = 4
Package Information:
Status `/builds/JuliaGPU/DiffEqTutorials.jl/tutorials/exercises/Project.toml`
[2169fc97-5a83-5252-b627-83903c6c433c] AlgebraicMultigrid 0.3.0
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.5.0
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.22.0
[aae7a2af-3d4f-5e19-a356-7da93b79d9d0] DiffEqFlux 1.16.0
[071ae1c0-96b5-11e9-1965-c90190d839ea] DiffEqGPU 1.3.0
[9fdde737-9c7f-55bf-ade8-46b3f136cc48] DiffEqOperators 4.10.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.14.0
[587475ba-b771-5e3f-ad9e-33799f191a9c] Flux 0.10.4
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.21.0
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 1.5.1
[47a9eef4-7e08-11e9-0b38-333d64bd3804] SparseDiffTools 1.9.0
[684fba80-ace3-11e9-3d08-3bc7ed6f96df] SparsityDetection 0.3.3
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 4.2.4